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Abstract 

We introduce a "loosely coherent" method for detection of contin- 
uous gravitational waves that bridges the gap between semi-coherent 
and purely coherent methods. Explicit control over accepted families 
of signals is used to increase sensitivity of power-based statistic while 
avoiding the high computational costs of conventional matched filters. 
Several examples as well as a prototype implementation are discussed. 

PACS numbers: 07.05.-t, 07.05.Fb, 04.80.Nn, 95.55.Ym 

1 Introduction 

The need for methods described in this paper arose during development of 
the PowerFlux search pLj for continuous gravitational wave signals. Even 
though aimed at a specific purpose of following up PowerFlux outliers, they 
have much wider applicability. To that end we will present a simplified 
description that omits some technicalities specific to searches for continuous 
gravitational waves. 

The PowerFlux algorithm [2] detects gravitational waves by computing 
power received from a particular direction at a certain frequency and spin- 
down. Similar approaches include Hough and StackSlide searches [U [31 HI |5l 
El d E] • Also, searches have been carried out with algorithms using substan- 
tially larger coherence lengths such as J-'-statistic P, [101 ITT] . 



The power-based methods are computationally efficient and allow all-sky 
blind searches to be performed with the sensitivity scaling as fourth root of 
the amount N of analyzed data. In contrast, coherent searches scale as A^~^/^ 
but become impractical for moderate values of A^. They also rely heavily on 
knowing the exact form of the expected signal - an assumption that we feel 
is overly bold when one is looking for a form of radiation for which no prior 
direct observation exists. 

There are searches that fill the space between these extremes. One way 
is to combine incoherently an output of multiple coherent searches. An- 
other approach is to perform a hierarchical search that follows up outliers 
with longer baseline coherent investigation. Both employ longer coherence 
baselines than power-based methods. 

Thus, in order to make a successful detection, one needs to overcome a 
"potential barrier" in computational costs that separates a blind search from 
an easy verification of a successful candidate. 

One reason for difficulties with current coherent methods is that they 
are optimized with a specific signal waveform in mind, and then the search is 
iterated over many signal templates. The templates often overlap [12j and, in 
fact, oversampling is routinely used to ensure that no signals are missed. This 
design is well warranted if sufficient computational power exists to exhaust 
the entire search space - but this is a situation current gravitational wave 
searches are not in. Furthermore, maximization alone is not necessarily the 
most optimal statistic [131 E] . 

We believe that an approach that combines attention to sensitivity and 
computational efficiency with more agile control over accepted waveforms is 
both more physically prudent and computationally accessible. To illustrate 
this, we present a loosely coherent method that is based on estimating power 
for a family of signal waveforms at once. 

2 Statement of the problem 

For the purposes of this paper we will assume that our entire dataset has 
been broken up into short portions each of which has been subjected to 
the Fourier transform, and we are looking for a signal of constant amplitude 
that would land into a single frequency bin {ak}^^i in k-th short Fourier 
transform with varying phases {(l)k}k=i- 

If the phases were known in advance we could compute the power of a 
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the high values of which would indicate the presence of the signal. There 
is a large body of literature that describes designing statistics with optimal 
signal-to-noise ratios (SNR), in particular [15] . 

In many cases a part of signal evolution (such as Doppler modulation 
induced by motion of the Earth) is known in advance. If we assume that this 
contribution has been factored out then the coherent power sum reduces to 
the case 0/+! — 0; = 0. 

The set T of all possible phases (modulo 27r) forms an A^- dimensional torus 
on which P is a smooth function. In practice, phases cannot be determined 
exactly ahead of time, but rather obey a set of constraints. Such family of 
signals would sweep a submanifold S* C T, possibly with boundary. 

Our goal is then to find a statistic that achieves high values when a 
signal from S is present and low values otherwise. One way to do that 
is take the maximum of P over {(pk} constrained to the submanifold S. 
Another approach is to view the unknown parameters as random, with the 
phases forming stochastic process, usually highly correlated. It is important 
to note that for either detection or establishment of upper limits we only 
need to know whether the signal is present, as the parameter estimation can 
be performed by partitioning S into subsets. 

We call this a loosely coherent approach, as instead of trying to find signals 
with a certain pre-determined set of phases, we are content with any signal 
that has phase evolution from S. The choice of the set S and the statistic P is 
then up to the designer of the search thus providing the necessary freedom to 
satisfy conflicting demands of efficiency in computation and signal recovery. 

Of course, any practical detection algorithm, even designed with full 
knowledge of expected signal, will respond to data with signals from a wider 
set of phases than physically expected. Tailoring the set S at the design 
stage, rather than simply characterizing it after implementation, allows finer 
control over which astrophysical signals one can detect and particulars of 
template placement. 



3 



3 Implementation of loosely coherent statis- 
tics 



3.1 Maximization 

The most straightforward way to construct a loosely coherent statistic is to 
maximize P over the set of possible phases S. This is a classical optimization 
problem with a quadratic objective that possesses several difficulties: 

First, we are trying to maximize a non- negative definite quadratic func- 
tion - thus our problem is inherently non-convexQ even for small portions of 
S. This precludes the use of well known optimization methods like gradient 
descent. 

Secondly, the dimension is very large, with small searches starting at 
= 1000. 

The third difficulty is more subtle and is due to the nature of interesting 
signal families S. These usually involve phases that evolve moderately fast 
with k and can wrap around numerous times. A typical example is a linear 
evolution produced by mismatch in frequency given by 

(j)k = A + Bk (2) 

with B on the order of 0.1. 

Because of the wrap around, a small uncertainty in (pk for some k can 
result in very large uncertainty in (pi for |Z — A;| 3> 1. In the limit N ^ oo 
the embedding of S into the torus T°° (considered with L°° norm in which 
it is not compact) stops being differentiable or continuous altogether. 

The properties of the map : 5* — )■ T as A^ approaches infinity are tightly 
connected with the scalability in the number of templates. To describe this 
connection we need some well-known tools from functional analysis. 

Let S" be a bounded (i.e. compact) finite dimensional manifold, possibly 
with boundary, with a metric ps- As mentioned before, we consider the torii 

with metric 

PJv({0fc}f=i, {■ipk}k=i) = supinf -ipk- 2vrm| (3) 

k m 

maximization problem maxx&s f{x) is called convex if the set of points 
{{x,y) : X G S and y < f{x)} is convex. In particular, for a differentiable /, this assures 
that the gradient descent method can not become stuck in a valley. 
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Let $Ar : S" — )■ be the family of embeddings describing phase evolution 
for successive SFTs. 

Our goal is to select templates in S such their image under $ ^ forms an 
e-net - any point in ^^{S) is within e of an image of some template. 

We distinguish three fundamentally different situations: 

• The map $00 : "S* — t- T°° is Lipschitz, i.e. it satisfies the following 
property: 

p^mxo), <^oo{xi)) < Lpsixo, xi) (4) 

Any continuously differentiable map is Lipschitz. In this case, we can 
cover $00(5') with any desired tolerance e by constructing a set of tem- 
plates in S which forms an e/L-net. A well-known fact from topology 
|16j is that it is possible to find coverings with template count scaling 
as e""^ where d is the Hausdorff dimension of S. 

Thus, we see that the template count does not depend on N and is 
proportional to £-'^^(■5') _ the best we could hope for. An example of 
such a map is given by 

^^{A,B) = {Asm{uk + B)}'^^, (5) 

where a; is a fixed parameter (such as Earth rotation frequency) and A 
and B are bounded search parameters. A physically relevant example 
is given by phase shifts from amplitude response of the detector. 

• The map $00 : S — ?■ T°° is known to be continuous, but not Lipschitz. 
In this case, we can still find a suitable template set for any desired tol- 
erance e, but the spacing of the templates in S will not depend linearly 
on e as it does in the Lipschitz case. We thus retain independence of 

but the number of required templates can grow faster than e-dim(5)_ 

A mathematical example of such a map is given by 

The required template count grows as e~^. We are not aware of any 
physically motivated search for continuous gravitational radiation that 
has parameters of this form. 
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• The map '■ S T°° is not continuous. While this can be due 
to trivial causes such as partial breaks in otherwise Lipschitz map, in 
general it would not be possible to find a finite template set to cover 
T°°. For the finite case $00 '■ S — )■ the template count will grow 
with N. 

An example of such a map is given by frequency evolution discussed 
above: 

= {Ak}l, (7) 
for which the required template count scales as A^e~^. 

One way to deal with these difficulties is to partition into small enough 
sets so that maximization can actually be carried out and combine the results 
afterwards. Further computational savings result from picking 5* described by 
only a few necessary parameters and overcoming their scaling properties with 
large computing power. The coherent searches for gravitational radiation 
such as P [ini [H] can be viewed as examples of this approach. 

3.2 Averaging 

Another way to bring computational costs under control is to replace P 
with a related function with a smaller Lipschitz constant. One can achieve 
this by averaging P over S or its subsets, which is equivalent to computing 
expectation value of P over some assumed distribution on S. This spreads the 
signal response over a larger area, but we only have to make the computation 
once for each subset. For ease of exposition we use the usual Lebesgue 
measure and average power rather than a more complicated statistic such as 
likelihood. 

In the most extreme case we just average away the phases (pk yielding the 
conventional semi-coherent method: 

N N 

EP = E^4a/e'('^'=-'^') = ^|a/p (8) 

k,l=l 1=1 

If the phases are truly random this statistic will perform better in the pres- 
ence of well behaved noise than computation of the maximum [T5] . 
A more conservative approach will limit phase evolution: 

S = {{M : \<Pi - <l>i+i\ < S} (9) 
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yielding the following statistic (computed using variables Si = (pi^i — (pi): 



EP=pp^ j ^..J ^PdS,...dS^^^=J^alail^^j (10) 

which interpolates between the fully coherent sum for 5 = and the semi- 
coherent case 6 = TT. The allowed spacing between frequency templates 
increases with 6, and in the limiting case — )■ oo is determined by the 
value of 6 /it in units of frequency bins. This has proven to be a good initial 
estimate of the spacing required by searches where N ^ 1/6. 

This method will lose some power if the true frequency of the signal 
at the time corresponding to coefficient is not a harmonic sampled by 
the Fourier transform. To avoid this, one can replace with more precise 
values estimated from the Dirichlet kernel. This effectively makes sure that 
the point with all phases belongs to - a condition we assume from now 
on. 

It is also possible to use the same approach to reduce the influence of 
periodic changes of underlying frequency, such as caused by mismatch in sky 
position and the resultant Doppler shifts. Assume 

S = : h = A + Bsm{cok + C)} (11) 

where A is some unknown (and irrelevant) phase, cu and C are known and 
fixed (such as from sidereal Doppler modulation) and B is allowed to vary, 
subject to \B\ < p. Then 

EP = ^ I^, Em=i alaie^(^'^~MdB = 

k,l=l ^k^l l5{sin{uik+C)-sin{uil+C)) ~ \^'^) 

EN ^, sin(2/3 sm{uj(k-l) /2) cos(C)) 
k,l=l^k^l 2/3sin(w(fc-/)/2)cos(C) 

As we have chosen a simple power sum P as a starting point, our averaged 
statistic will always have the form 

N 

EP=J2 ('l^'iKki (13) 

k,l=l 

for some kernel K^i and is thus similar to cross-correlation search [17]. As 
we will see later the efficient computation of the sum for small 5 is best done 
in a manner different from the cross-correlation statistic. 
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3.3 Loosely coherent searches as a filtering problem 



The statistic EP can be rewritten as a scalar product of the vector of input 
data a with the image of a under the operator K' which square K' K' is 
given by the kernel i^fe^: 

N 

EP = ^ alaiKki = a^K'^R'a (14) 

k,l=l 

From this point of view K' acts as a filter rejecting signals outside the ex- 
pected set, after which we take the usual semi-coherent sum. 

For example, K' can be chosen as a low-pass filter given by a sine or 
Lanczos kernel. This would admit signals with phases varying slower than 
the filter cutoff frequency. 

For a practical implementation the main point of concern is the ability 
of the statistic to tolerate frequency mismatch, as it directly impacts the 
number of templates. For this purpose the low pass filters are optimum, 
tolerating mismatch values up to a cutoff frequency and rejecting signals 
with faster varying phases. 

A more sophisticated approach is to assume a distribution on the set of 
allowed phases and then treat our signal as a highly correlated stochastic 
process. Since the data analysis is typically carried out after the data collec- 
tion is complete, one is not restricted to causal filters alone and, in the case 
of stationary noise and limited phase evolution, we obtain a low pass filter 
as a solution. 

The loosely coherent statistic based on a sine filter is optimal in the 
following idealized situation: suppose our data {a^} consists of a sum of 
stationary mean zero Gaussian noise of known variation (which is typically 
easy to estimate from data known not to contain any signals) and unknown 
band-limited signal of limited power, with no additional information on the 
signal form or phase evolution. A Fourier transform will separate our data 
into high frequency area where there is no signal and which can be safely 
discarded and low frequency area which phase information is irrelevant due 
to the signal having an arbitrary spectral shape. 

We are thus left with a problem of deciding whether our low-frequency 
data is consistent with Gaussian noise alone or there is an arbitrary additive 
signal present. 

Both the limited power condition and the structure of Gaussian noise 
are symmetric under unitary transformations. Thus, if no other restrictions 
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are present, the only meaningful information is the power contained in the 
low-frequency data. 

While this fairly standard argument bridges both frequentist and Baycsian 
approaches, it does have a number of limitations. The most severe is that 
the symmetry is lost in case of non-stationary noise. Additionally, a family 
of physical signals can be expected to have a spectrum more interesting than 
a plain flat-top. 



4 Practicalities of the gravitational wave searches 



We will now qualify the phase shift evolution that one expects to encounter 
in current searches. 

At the moment, the searches analyze data from 50 Hz through 1500 Hz, 
accounting for spindowns as large as — 10~^ Hz/s. The analysis is done using 
short Fourier transforms (SFTs) of ^ 1800 s length, which have 50% overlap 
in some searches, no overlap in others and often have gaps. For this paper 
we will assume that the time interval At between and a^+i is 1800 s. 

We will assume that {ak} have already been adjusted so that the template 
O with all (pk — is in S. 

There are several sources of non-trivial phase shifts, which we will describe 
in terms of maximum expected difference 6 between nearby phases: 

• Frequency mismatch - a template possessing frequency different from 
O by A/ will experience a linear phase evolution of 



• Sky position mismatch - a mismatch in sky position will produce a 
slightly different Doppler shift. On short time scales this is dominated 
by Earth rotation (with velocity ~ lxlO~^c) and is periodic in time 
and linear in sampled frequency: 



where Ar is the maximum expected mismatch in radians, with practical 
values usually less than 0.01. 



5 = 27rA/Ai 



(15) 



5 = 27rlO-VAiAr 



(16) 
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Spindown mismatch - a spindown different from O by A/ will produce 
a linear evolution of the frequency and, thus, a quadratic change in 
phase: 

5 = 27rA/TAt (17) 

Here T shows maximum variation of time variable with respect to ref- 
erence time. If the reference time is positioned at the center of the run, 
then T is half the time base. 

Source frequency evolution - the source signal can be modulated by 
a nearby orbiting object. Assuming circular orbit with radius r (ex- 
pressed in astronomical units) and using p = m/M for the ratio of 
object mass m to the star mass M (both expressed in units of solar 
mass) the angular frequency of the modulation is: 



and the maximum Doppler shift from the central body is 



c c \ r \ r 

The worst case change in phase induced by this motion over time At 
and assuming radiating frequency / is: 

S = /HH.At « 2.2x10- Hz.^ ^i^^pyiT^ (20) 

c 1000 Hz 1800 s(r/l AU)2^^ ^ ^ 

The curiously small size of 6 is due largely to the small value of product 
uAt. For a search that assumes a specific phase evolution over a long 
time interval this would be much larger. The loosely coherent search 
is not completely immune from this effect - it will lose power when 
enough phase accumulates during integration for the signal to escape 
into nearby frequency bin. This suggests that searches looking for 
more extreme systems should use coarser frequency bins, smaller At 
and tighter 6. 

Table [T] shows the expected phase shift for conditions commonly encoun- 
tered in present day searches. 
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Table 1: Maximum phase change in degrees between 1800 second spaced 
SFTs. 



Phase shift cause 100 Hz 500 Hz 1000 Hz 2000 Hz 

Frequency mismatch of 36 36 36 36 

A/ = 0.1/At 

Sky position mismatch of 1.1 6 11 23 

Ar = 1° 

Spindown mismatch of 20 20 20 20 

A/ = 10-12 H^/s for T = 1 y 

Source modulation for 0.1 0.6 1 2 

p = 1 and r = 0.1 AU 



5 Efficient computation of loosely coherent 
sums 

We will now turn to efficient computation of the loosely coherent statistic. 
Given reduced sensitivity to perturbations in search parameters compared 
with purely coherent methods and corresponding reduction in the number of 



templates, the quadratic cost of computing the sum (13) is not completely 
unreasonable. 

Noticing that the kernel Ki^ is a positive symmetric matrix, one expects 
to do better by finding eigenvectors and eigenvalues of K and discarding 
eigenvectors with small eigenvalues. This will make the computational cost 
bilinear in and the number of remaining eigenvectors. 

Let us consider, as an example, the case of limited phase evolution — 



>k+l 



< 5 with the previously computed kernel 



A-,.„ = ('E^l^V'^L e-"K-l (21) 



where we introduced a = — log 

When 5 = we are dealing with a fully coherent case and the kernel has 
only one eigenvector with non-zero eigenvalue, while for 5 = oo we have the 
semi-coherent case and K is the identity matrix for which we have to use 
the entire basis. It seems reasonable to expect that for small 5 we will have 
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a few-eigenvector situation, while for large S we will have something similar 
to a semi-coherent sum, where it makes sense not to truncate by eigenvalue 
but rather cut side diagonals of K that are small. 

It turns out that the set of "small" 6 values is quite large. To see why 
this is so, first examine the plot of a versus phase mismatch 6 on figure [T] 
Even for a phase mismatch as much as 45° the value of a is relatively small 
at ^ 0.1. 



s 




phase mismatch (degrees) 



Figure 1: Dependence of a on phase mismatch 6 



Secondly, consider the continuous version of our kernel: 

/oo 
e-"l^-^l/(t;)rft; (22) 
'OO 

The operator K is given by a convolution of f{v) with e~"'"'. As is well- 
known, Fourier transform will convert convolution into multiplication. Thus 
the spectrum of the convolution operator is given as Fourier transform of 
its kernel. The functions e*^" can be considered as eigenvectors of K in 
appropriate functional space (e.g. C°°): 

^-a\u-v\^iXv^^ = — ^^^e'^" (23) 



u 
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The eigenvalues have the famihar Lorentzian form ^^^^ with quadratic de- 
cay. In hindsight, this is not surprising as the condition \(f)k — <Pk+i\ < S 
is similar to the requirement that the signals we are looking for are band 
limited. 



Table 2: Number of eigenvectors required to compute K with 1% accuracy. 



Phase shift 6 N = 1000 N = 5000 N = 10000 



1° 


2 


3 


4 


5° 


7 


22 


42 


10° 


18 


81 


162 


20° 


66 


324 


647 


30° 


148 


737 


1474 


45° 


351 


1751 


3502 



Table 3: Number of eigenvectors required to compute K with 5% accuracy. 



Phase shift 6 N = 1000 = 5000 = 10000 



1° 


1 


2 


2 


5° 


3 


10 


18 


10° 


8 


36 


71 


20° 


29 


142 


283 


30° 


65 


321 


641 


45° 


147 


736 


1471 



Tables [2] and |3] show the number of eigenvectors needed to approximate 
K given by formula 21 for various numbers of equally spaced SFTs and 
some typical values of S. The approximation is done using operator norm 
which is equivalent to counting the number of eigenvalues that are at least 
1% for table |2] (5% for table [3]) of the largest eigenvalue of K. While the 
fraction of eigenvectors does rise linearly with and thus the computational 
requirements are still quadratic, said fraction is a rather small number for 
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S < 10° and in practical implementations (especially on processors with 
vector arithmetic) the scaling will be close to linear. 

For larger values of 6 one might wish to go with a different algorithm. In 
particular, it makes sense to consider decompositions using non-orthogonal 
vectors, the simplest of which is obtained by truncation of side diagonals. 



E3 o 
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Figure 2: Eigenvectors of the kernel 21 corresponding to first four largest 
eigenvalues for = 1000 and 5 = 5° 



As we mentioned before, in the continuous case the eigenvectors are simple 



sme waves e 



The discrete case is nearly sinusoidal. The eigenvectors of 
the kernel 21 corresponding to first four largest eigenvalues for = 1000 
and 5 = 5° are shown on the figure |2} The eigenvector corresponding to 
the largest eigenvalue is not constant and can be regarded as a window one 
applies to the data in order to make the usual coherent sum respond to signals 
from S. The eigenvector decomposition was done numerically using R [TH] . 

This idea can be exploited to speed up eigenvector decomposition, by 
analytically transitioning into the basis of pure sine waves and then discard- 
ing entries of K from higher order modes. The remaining matrix of smaller 
dimension can then be diagonalized with conventional numerical techniques. 
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It is interesting to consider the case of very short Fourier transforms of 
a few seconds in length and correspondingly small At. The phase shifts 
from most sources (except for frequency mismatch) will be small as well, 
and computation of K can be performed by taking a Fourier transform of 
the input data and then summing up power in low frequency harmonics 
weighted by eigenvalues of K. 

This has close relation to the resampling technique [19]. 

The resampling implementation of J-'-statistic operates by heterodyning 
30 minute SFTs to a desired frequency, inverting the Fourier transform to 
obtain a time series which is stitched together and then band-limited and 
downsampled. The resulting time series is converted into detector frame 
which allows efficient computation of J-'-statistic using Fourier transform. 

Another way to obtain the same time series is to start with shorter SFTs 
which frequency bins are large enough to accommodate Doppler shift. A 
time series of frequency bins of these short SFTs is then just another way of 
heterodyning our input data with the advantage of bypassing the need for 
inverse Fourier transform. If the frequency band that is being searched is 
significantly smaller than the size of initial frequency bins the time series can 
be band-limited and downsampled just as done in [19]. 

The conversion of heterodyned time-series into detector frame consists of 
two parts: removal of the phase shift from signal evolution due to intrinsic 
effects or Earth motion, which is also done by loosely coherent method, and 
interpolation in order to obtain evenly spaced time series suitable for fast 
Fourier transform algorithm. 

The computation of J-'-statistic involves summing three terms quadratic in 
the elements of our time series with coefficients that depend on time position 
of the source and the detector but not the amplitude or polarization of the 
expected signal. This can be viewed as computation of a specific kernel K 
which rank is at most 2. If we take the interpolation algorithm into account 
the rank will increase but will still be much smaller than kernel dimension. 

The same approach can be used to compute loosely coherent statistic 
where we might need to use additional terms to accommodate kernels with 
larger rank. In return, the statistic can be made more tolerant of mismatch 
in source parameters, such as sky location. 
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6 Sensitivity estimates 



It must be said that the sensitivity of a given method is best judged from 
a search made on real data, as computational efficiency and practicalities of 
detector artifacts in the input data have often a much stronger impact than 
an extra few percent gained by fine-tuning the algorithm with analytical 
considerations that assume Gaussian noise. 

Nevertheless, it is useful to have an idea of what to expect in the perfect 
situation as a starting point for practical applications. We will concentrate on 
the case of perfectly coherent signal and how the performance varies between 
the extremes of coherent and semi-coherent power sums. 

The standard methods of filtering theory can be employed to obtain a 
rough estimate. As we mentioned before, the phase evolution condition — 
(pk+i I < ^ is closely related to the condition that our signals are band limited. 
In this case, the rejection of noise outside the acceptance band results in 
improvement in the signal-to-noise ratio compared to the usual semi-coherent 
case which is sensitive to all signals within the frequency bin of the original 
SFTs. 

The acceptance band is narrowed down by a factor inversely proportional 
to the number of SFTs it takes for the phase to make a full turn (not to 
exceed, of course, the total number of SFTs available). Thus, given a fixed 
number of SFTs, we expect the improvement in the signal-to-noise ratio to 
scale as 1/ \/6 tempered by the non-linear effects of our statistic. 

This is illustrated on figure [3] that shows results of simulation evaluating 
signal-to-noise ratio gain for limited phase evolution statistic as we decrease 
S for a coherent signal. The simulation was performed using = 1000 
SFTs which were composed of Gaussian noise C,i with standard deviation 1 
and a constant signal with amplitude h = 0.7 which results in the average 
signal-to-noise ratio of ~ 8 for a semi-coherent search. 

The statistic was computed according to the formula 

N 

P{h) = Y,K,,{i, + h){l, + h) (24) 
where the kernel Kij was either an identity matrix for semi-coherent case, a 



matrix with 1 in all cells for the coherent case or given by the formula |2T] for 
the loosely coherent case. 
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The signal-to-noise ratio in this simulation was defined as the value of the 
statistic minus the average value obtained on noise alone and divided by the 
standard deviation of values produced by pure noise: 

mean(P(/i)) — mean(P(0)) 



SNR 



sd(P(0)) 



(25) 



Here mean and standard deviation were taken over 1000 independent real- 
izations of noise. 

All of the statistic values are described by a weighted x-squared distribu- 
tion which depends on 6. For large 6, however, it is close to a Gaussian dis- 
tribution as well due to the central limit theorem. To illustrate the change in 
the distribution of our statistic we show 10% and 90% quantiles of the signal- 
to-noise ratios obtained as well as the mean. The vertical axis is logarithmic, 
so the spread in signal-to-noise ratios increases as 6 becomes smaller. 

q10 o mean o q90 o 



300 
200 



cr 50 

z 

(n 

30 
20 



10 
7 

6 -I 




coherent pi/100 pi/50 pi/25 pi/10 pi/4 pi/2 semi-coherent 
5 



Figure 3: Dependence of signal-to-noise ratio on phase mismatch S. The 
upper, central and lower curves show 90% quantile, mean and 10% quantile 
of multiple simulation runs. 



The fiattening out of the curve for small 5 is due to different scaling 
regimes near the extremes of coherent and semi-coherent statistics. This can 
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be illustrated by considering a semi- coherent statistic that operates on 
SFTs which are coherently combined in stretches of k SFTs each and the 
results are combined incoherently. Then the scaling law for the signal-to- 
noise ratio is: 

SNR ~ k^/N/k (26) 
Now suppose that k = aN is a certain fraction of N. Then the scaling is 

SNR ~ (27) 

As our statistic is power based the sensivity will scale as l/{^/a\/N). The 
fourth root in a has a really slow growth. For example, for a = 0.1 it is only 
0.56 - so for less than a factor of 2 loss in sensitivity the coherence length 
can be dropped by a factor of 10. 



7 Prototype implementation 

An initial implementation of the loosely coherent statistic was done within 
the framework of the PowerFlux [2] program. This implementation provided 
practical experience with a loosely coherent search and addressed the problem 
of following up outliers from the all-sky PowerFlux search over LIGO's fifth 
science run. 

As the underlying code base was not designed with the loosely coherent 
search in mind, the code has a number of inefficiencies. In particular, the 
double sum in the statistic EP was computed by brute force. Nevertheless, 
the speed was sufficient to quickly carry out searches in disks of 0.03 radians 
radius on the sky over 20002 SFTs split evenly between HI and LI detectors. 
The powers from individual detectors were combined incoherently to make 
the comparison to semi-coherent code more fair. The nearby SFTs were 
separated by 30 min. In practical data, the SFTs are usually 50% overlapped, 
but there are can also be gaps in the data. The 30 min constant was chosen 
reasonable worst case. 

While the analysis of actual interferometer data is still underway, we can 
report on results of simulations using Gaussian data. For these simulations 
we used a Lanczos kernel with parameter 3: 

sin(<5(t2-ti)/30 min)sin(i(5(t2-ti)/30 min) , <5|t2-ii| 

^ ' M when > 37r 

30 mm — 
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This kernel naturally vanishes for widely separated SFTs which makes this 
a variant of cross-correlation search, albeit with particularly large number of 
off-diagonal entries, which is further increased by the 50% overlap of nearby 
SFTs that is usually employed by PowerFlux. We explored values of 6 as 
small as 7r/5 which involves summing up to 59 diagonals when working with 
overlapped SFTs. For these values of 6 the required computational time 
scales as square of observation time (for time bases several months and larger) 
and cube of covered frequency range. 



semicoherent search + 
pi/2 loosely coherent search o 
injected value 




0.0 0.1 0.2 0.3 0.4 0.5 
Frequency mismatch as fraction of bin size 



Figure 4: Upper limit versus frequency bin mismatch. The horizontal line 
marks the strain of the software injections. The upper curve shows upper lim- 
its from a semi-coherent search which are consistently above injected value. 
The upper limits from loosely coherent search follow semi-coherent search 
before the limit of phase tolerance is reached and decline sharply afterwards. 

Figures |4] and |5] show results of Monte-Carlo injection run assuming a 
static source location (right ascension 2.0, declination 1.0, spindown 0) and a 
linearly polarized signal. This choice was made to increase readability of the 
plots as all-sky injections with arbitrary polarizations inject different amount 
of power in the interferometer making the curves wider. The injections were 
made into Gaussian data that was filtered to simulate Hann windowed short 
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Strain of Injected signal 



Figure 5: Signal to noise ratio versus strain of Monte-Carlo injections. The 
values were capped at 15 in order to expose the more interesting low SNR 
region. 



20 



Fourier transforms (SFTs). The assumed frequency range varied from 400 
to 410 Hz and SET frequency bin size was 1/1800 Hz. 

The 95% confidence level upper limits are produced by PowerFlux code 
for a set of 501 frequency bins given a particular direction on the sky and a 
spindown value. The results are then maximized over a set of polarizations 
and small area on the sky around the injection point. This follows the analysis 
method used in [1] and [20] . 

Both semi-coherent (power only) and loosely coherent algorithms pro- 
ceed by sampling discrete range of frequencies with configurable spacing in 
fractions of SFT bin size. Figure |4] compares how the mismatch between the 
actual injected frequency and the sampled frequency affects upper limits pro- 
duced by semi-coherent and loosely coherent codes. The frequency spacing 
was set at 1 SFT bin and the injected strain value was fixed to 1.8x10"^^. 
We see that a loosely coherent search with 5 = tt/2 has an initial fiat response 
for small mismatch in frequency which is followed by rapid decay to values 
below injected strain. In contrast, the semi-coherent search shows only minor 
reduction in the upper limit which is fully compensated by built-in correction 
factor. 

Figure [5] compares the signal to noise ratios (SNRs) of semi-coherent and 
loose-coherent methods. The frequency spacing of the loosely coherent search 
was reduced to l/8th of the SFT bin which insures correct reconstruction of 
the upper limit for the entire range of weak and strong signals. Because of 
the larger number of templates, the SNR achieved on pure noise is higher 
for the loosely coherent search than that of the semi-coherent search. For 
signals above noise the loosely coherent search produces signal-to-noise ratios 
on average 50% larger than semi-coherent one. 

8 Summary 

We have discussed the problem of detecting a family of signals S from the 
point of view of computational efficiency and presented a method of creating 
a statistic that is sensitive to the entire family S or its subset. Two simple 
examples were considered which showed close ties to well-known methods of 
matched filtering, cross-correlation and semi-coherent sums. 
There are several directions of further study: 

• The prototype large 5 implementation shows feasibility of the overall 
method, but does not provide information on the overall computational 
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efficiency. We plan to develop a dedicated small S code to be used in 

targeted searches that cover small sky area (such as galactic center 
or globular clusters). This should provide experience with scalability 
properties of the loosely coherent method. 

• The average of P was used to make the maximization computationally 
tractable. In fact, for small the maximization can be carried out 
directly. It is worthwhile to investigate the possibility of combining the 
two techniques. 

• For the case of the set S given by conditions \4>k — 4>k-i-i\ < S and 
assuming small 5 the maximization over P can be carried out assuming 
(pk+i — 0fe i This converts the problem into the discrete domain and 
makes it amenable to binary optimization methods which have seen 
much progress in recent years. A particularly interesting observation 
is that for a noise dominated signal the function to be optimized has 
random coefficients, so an optimization method that works only on a 
certain proportion of objective functions can yield useful results. 
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